PERIODIC OSCILLATIONS OF BLOOD CELL POPULATIONS IN 
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Abstract. We develop some techniques to prove analytically the existence and stability of long 
period oscillations of stem cell populations in the case of periodic chronic myelogenous leukemia. 
Such a periodic oscillation pea can be analytically constructed when the hill coefficient involved in 
the nonlinear feedback is infinite, and we show it is possible to obtain a contractive returning map 
(for the semiflow defined by the modeling functional differential equation) in a closed and convex 
cone containing poo when the hill coefficient is large, and the fixed point of such a contractive map 
gives the long period oscillation previously observed both numerically and experimentally. 
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1. Introduction. "How do 'short' cell cycles give rise to 'long' period oscilla- 
tions?" This question has arisen from the observation of blood cell population oscil- 
lations in the case of periodic myelogenous leukemia (PCML) jfT,, a blood disease to 
be discussed in some details below. Indeed, it has long been observed in the bone 
marrow that there is an enormous difference between the relatively short cell cycle 
duration which ranges between 1 to 4 days \T2\, > EE] an d the long period oscil- 
lations in PCML (between 40 to 80 days) [J]. The link between these relatively short 
cycle durations and the long periods of peripheral cell oscillations is unclear, to the 
best of our knowledge, has neither been biologically explained nor understood. An 
attempt to answer this question has been made by Pujo-Menjouet and Mackey in 
and Pujo-Menjouet et al. in [22], where they investigated the role of each parameter 
of the mathematical model involved in the cell cycle and the influence of each param- 
eter on the long period and the amplitude of the peripheral cell oscillations. They 
showed qualitatively that the cell cycle regulation parameters have major influence on 
the oscillation amplitude while the oscillation period is correlated with the cell death 
and differentiation parameters, and they obtained these results in the particular case 
where the hill coefficient involved in the model formulation is infinite. Our objective 
here is to prove analytically that the similar conclusions and results remain true in 
the more biologically realistic case where the hill coefficient is finite. 

More specifically, from the previous studies, it is known that the evolution of 
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the cells in resting phase involves the Hill function in both the term representing the 
instantaneous loss of proliferating cells to cell division and to differentiation and the 
term representing the delayed production of proliferating stem cells. A key parameter 
in the Hill function is the integer n which is usually large, and this Hill function 
reduces to the Heaviside step function when n = oo. As will be shown, the underlying 
system with n — oo becomes a piecewise linear scalar delay differential equation 
that, after non-trivial but straightforward calculations, has a periodic solution of 
large periods and amplitudes with very strong stability and attractivity properties. 
The main purpose of this paper is to construct a convex closed cone containing the 
aforementioned periodic solution (when n = oo) and a contractive returning map 
defined on this cone such that a fixed point of such a returning map gives a stable 
periodic solution of the model equation for the cells in resting phase when n is large. 
This method was first developed in Walther (2H], [HOI f° r a scalar delay differential 
equation with constant linear instantaneous friction and a negative delayed feedback, 
and was later extended to state-dependent delay differential equations [31 ,[32 and 
to delay differential systems (H2], ESI This method was further enhanced recently 
in j2l] by incorporating some ideas from classical asymptotic analysis and matching 
method. Applications of this method to the model for cells in the resting phase seem 
to be highly non-trivial since both the instantaneous loss of proliferating cells and the 
delayed production of proliferating stem cells involve the nonlinearity and there is no 
analytic formula for the periodic solution in the limiting case (n = oo). 

We should emphasize that periodic hematological diseases have attracted a signif- 
icant amount of modelling attention in various domain such as periodic auto-immune 
hemolytic anemia ^Hj an d cyclical thrombocytopenia |2S], (23- It has been ob- 
served that the periodic hematological diseases of this type involve periodicity between 
two and four times the bone marrow production delay. This observation has a clear 
explanation within a modelling context. Some other hematological diseases such like 
cyclical neutropenia (0,0) |10j . |14j . jloj . \17\ ) and chronic myelogenous leukemia (3 
involve more than one blood cell type (i.e. white cells, red blood cells and platelets). 
It is believed that the oscillations in these diseases originate in the pluripotential stem 
cell compartment and have very long period durations (of order of weeks to months) 
in general. In the particular case of the periodic chronic myelogenous leukemia, the 
period can range from 40 to 80 days and two lines of evidence appear to prove that 
the oscillations are due to a destabilization of the stem cell population based in the 
bone marrow. The first evidence is due to the gene mutation in the Philadelphia 
chromosome and responsible for the disease . The mutated cells have been observed 
in all the blood cell lineages @], El- [HI, |1 1| • |28j . The second line of evidence is given 
in (3 where the authors collected clinical data from the literature, and proved that 
white blood cells, erythrocytes and platelets oscillate with the same period. 

Periodic chronic myelogenous leukemia (or PCML) takes its name from the clinical 
character and the type of leukemia it describes. Leukemia is a malignant disease 
characterized by uncontrolled proliferation of immature and abnormal white blood 
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cells in the bone marrow, the blood, the spleen and the liver. Its character can be 
chronic (in the early stage of myelogenous leukemia) or acute (in the late stage) . The 
type of cells involved are myeloid, lymphoid or monocytic depending on the damaged 
branch of blood cell production. The stem cell model by Pujo-Menjouet and Mackey in 
[23 and Pujo-Menjouet et al. in (221 (also called Go model due to the consideration of 
the Go resting phase in the cell cycle) is developed in order to describe the mechanism 
of the disease under consideration 0, El) EH- 

The remaining part of this paper is organized as follows. In Section (21 we present 
the model in detail. In Section Owe recall some previous results from Pujo-Menjouet 
et al. in [22] in the case where the Hill coefficient n is infinite. Then, we introduce a 
more general result on the perturbed delay equation given in sectional and we present 
our main results in section including the full asymptotic expansion for the periodic 
solutions. 

2. Description of the model. The Go model, whose early features are due to 
Lajtha El and Burns and Tannock 0, is derived from an age structured coupled 
system of two partial differential equations, coupled with some boundary and initial 
conditions [22], El, El> an d j20| . Using the method of integration along charac- 
teristics :83f these equations can be transformed into a pair of non-linear first-order 
differential delay equations El' El' El- The model consists of a proliferating phase 
where the cell population is denoted by P(t) at time t, and a Go resting phase, with a 
population of cell N(t). In the proliferating phase, cells are committed to undergo cell 
division a constant time r after their entry. Note that the choice of r as a constant 
is to simplify the problem, though some models with a non constant value of r exist 
PQ, P). The loss rate 7 in the proliferating phase is due to apoptosis (programmed 
cell death) . At the point of cytokinesis (cell division) , a cell divides into two daughter 
cells which enter the resting phase. In this phase, cells can not divide but they have 
the choice of between three different fates. They may have one of three possible fates: 
differentiate at a constant rate 5, reenter the proliferating phase at a rate 0, or simply 
remain in Go- Note that the reentering rate (3 will be a nonlinear term in our equation 
and the focus of our study (see Figure 112. lfl for an schematic illustration of the cell 
cycle). 

The model, described by a coupled non-linear first order delay equations, takes 
the following form 

^ = - 7 P(t) + (3(N)N - e-^(3(N T )N T , (2.1) 

and 

dN(t) 



dt 



-\P(N) +S}N + 2e~^ T p(N T )N T , (2.2) 



where N T = N(t — t). The resting to proliferative phase feedback rate j3 is taken to 
be a Hill function of the form 



4 



M. C. MACKEY, C. OU, L. PUJO-MENJOUET and J. WU 




Fig. 2.1. A schematic representation of the Go stem cell model. Proliferating phase cells (P) 
include those cells in G\, S (DNA synthesis), Gi, and M (mitosis) while the resting phase (N) cells 
are in the Go phase. 8 is the rate of differentiation into all the committed stem cell populations, and 
7 represents a loss of proliferating phase cells due to apoptosis. f3 is the rate of cell reentry from Go 
into the proliferative phase, and r is the duration of the proliferative phase. See \14\ , ]l f>f - f J *?j f for 
further details. 



In equation IJ2.2I) , the first term represents the loss of proliferating cells to cell division 
(J3(N)N) and to differentiation (SN). The second term represents the production of 
proliferating stem cells, with the factor 2 accounting for the amplifying effect of cell 
division while e _7r accounts for the attenuation due to apoptosis. Note that, we only 
need to understand the dynamics of the Go phase resting population (governed by 
equation l|2.2|l ) since the proliferating phase dynamics (governed by equation (j23J) 
are driven by the dynamics of the resting cells. 

By introducing the dimensionless variable x — N/9, we can rewrite equation 112.211 

as 

dx 



where 



\/3(x) + 5}x + k0(x T )x T , (2.3) 



and k — 2e 7T '. The steady state x* of equation 112.31 ) are given by the solution of 
dx/dt = 0. Then we have x* = , or 

k - 1 x 1/n 



= I /3 — ^— - I J ■ (2.5) 
Here we require 

T< -7 ln ^AT' 

k — 1 

so that /? — r - > 1 m ESI- 

o 

Note that when n — > oo, x* — > 1 in 112.5 1 and f3{x) tends to a piecewise constant 
function (the Heaviside step function). 
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A solution of Equation 112. 3|) is a continuous function x : [—r, +00) — > R + obeying 
^2.311 for all t > 0. The continuous function tp : [— r, 0) — > R + , </?(/;) = a;(t) for all 
i £ [— t, 0], is called the initial condition for x. Using the method of steps, it is easy 
to prove that for every tp £ C([— t, 0]), where C([— r, 0]) is the space of continuous 
functions on [— r, 0], there is a unique solution of equation l12..'ill subject to the initial 
condition ip. 

3. Periodic solutions: limiting nonlinearity. In this section we study the 
dynamics of equation l|2.3|) when (3(x) is reduced to the step function 

Pi?) - ( 2' x > V 

\ A), a; < 1. 

As in the paper by Pujo-Menjouet et al. [22] . we introduce two constants by 

a = /3 + £, r = 2/3oe" 7r = kf3 . 
Inserting the above step function f3(x) into equation \2.-\l . we have 



dx 
~dt 



—ax, < a; < 1 < x Tl 

—ax + Tx Tl < x,x T < 1, 
— 5x + rx T , < x T < 1 < x, 



(3.1) 



where x T — x(t — t). 

For the equation i'Uii . we choose the initial function t^(t) > 1 + ?/ for t £ [— r, 0) 
and f(0) — I + 77 where 77 is a small positive constant chosen later. We should remark 
that if we choose ip(t) < 1 + 77 for £ £ [— r, 0), the results and the techniques to be 
obtained and developed are similar. By the continuity of the solution x, we have from 
equation H3.1fl the existence of a ti such that x(t) and x(t — r) are greater than 1 for 
t £ [0, ti) and cc(ti) = 1. The solution then satisfies 

dx 

— = -dx, for * £ [0,t{\. (3.2) 

Thus solving the above equation, we can have x(t) = ip(0) e - 5t = (1 + T]) e - St . It 
follows that 

lny(O) _ ln(l + g) 
tx - ^— - s . (3.3) 

In the next interval of time, defined by + r], the dynamics are given by 



(3.4) 



The solution is then given by x(t) — e -0 ^^* 1 ) for t £ [ti, 4i+t] and satisfies cc(tx+r) = 
e _QT which is independent of the initial function <p{t). In other words, the solution 
destroys all memory of the initial data. 
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The solution in the next interval will be such that x, x T < 1. In order that 
equation (|3.1|l has periodic solutions, we should impose an extra condition on T and 
a so that 

-ax + Tx T >0. (3.5) 

Otherwise, if —ax + Tx T < 0, then the solution may tend to zero as t approaches 
infinity and thus we cannot expect periodic solution. In particular, if 

—ax + Tx T w 0, 

then the solution may stay below the line of x — 1 so long that the analysis becomes 
very complicated. This is also undesirable biologically since the period will be ex- 
tremely long. Note that for t e [h + r, t\ + 2r] , x(t - r) = e -a(*-*i-^)_ Then, from 

(jS3J) we have — = —ax + Tx T = —ax + re - " 1 "^ - * 1-7 ") which gives 



(t) = e -«(*-*i-T)( e -ar + r(t _ h _ r ))_ (3 g) 



X 

For the sake of simplicity, we impose an extra condition on T: 

r>max{-(e ttr -e- QT ),ae QT }, (3.7) 

T 

so that i|3,5JI holds when x < 1 (due to x(t — r) > e _QT ), and also x(t\ + 2t) = 
e- aT {e- aT + Tt) > 1 by 1 

Since x(t) is increasing in t G [t\ + r, ti + 2r], there exists a unique point £2 £ 
(ii + t, ii + 2t) so that x(i 2 ) = 1- Assume t 2 — ti + t + u, u e (0, r). Then from 113.611 
we have 

Tii. (3.8) 



The above equation 113. 8JI is a transcendental equation and cannot be solved explicitly. 

But we can expand e au 
approximated equation 



(au) 2 

But we can expand e au by Taylor's series, that is 1 + awH — and solve u by the 



2 

au + Y u2 ~e- aT + Tu. 



Next for t 6 [t2,ti + 2r], the dynamics are 



which gives 



^ = -Sx + Tx T = -Sx + re -«(*-^-), 

at 



s(t) = e- 5 ^*- t2 ){l - ile Q( * 1+T ^ 5t2 (e"** - e-* t2 )}. (3.9) 
Po 



1 Note that this condition allows us to get the shortest period length for the solution. In order to 
get longer periods, we should assume other conditions on V such that x(t\ + 2r) < 1, thus the slope 
of the increasing part of the solution would be less steep. 
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Finally for te[h + 2r, f 2 +t], 

— = —5x + Tx T 
at 



that is 



where 



-Sx + Ye-°'^- tl - 2T \e- ar + T(t -t x - 2r)), 



x {t) = e -*(*-*i-a^) [ X ( tl + 2r) + T - j(h + 2r))] 



= f e- aT + T(t - tl - 2r) - ^(t-t^) _ 

(o — a) \ o — a J 

We now claim that 

x{t) > 1, t € (t a ,i a +r]. (3.10) 

Indeed, at the point t2, x{t%) = 1, a;(<2 — t) > e~ QT . By l|3.7ll we have 

a/(*a) > Sx + Tx T > 0. 

Suppose, y way of contradiction, that there exists a point h G (i2,ta + t) such that 
a;(/i) = 1, < and x(t) > 1 for i G (t2, ft-). Then using equation ( 13 . If) we have 

by 

a/(/i) = -S + Yx{h -t)>-5+ Te~ aT > 0. 

This is a contradiction and our claim is true. 

After the time t% + r, both x\ and x are greater than 1, and the solution satisfies 

x' = -Sx(t) (3.11) 

and hence is decreasing. Therefore, there exists a point, say t = d so that x(d) = 1. 
Now we can use l|3.1()H and l|3.11H to choose a small positive constant rj such that at 
some point T x < d 

x(T x ) = 1 + 7?, x(T x + s) > 1 + r,, s e [-r, 0). (3.12) 

Actually, this T x is exactly the period of the solution x(t). Summarizing the above 
analysis, we have the following result: 

Theorem 3.1. Assume that x is the solution of { H. 1)) subject to the initial con- 
dition <j> > 1 + 7] where n is a small positive constant defined in IS. imp . Suppose that 
r satisfies ( \H. H\) . Then x is a periodic solution. 
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4.1. Perturbed delay equation. With the detailed analysis of the Go model 
when the Hill function reduces to the Heaviside step function, we can now consider 
the general nonlinearity from the viewpoint of regular perturbation. More precisely, 
we consider the perturbed problem 

^ = -\P(y)+6}y + kp(y T )y T , (4.1) 

i.e., we return to the original problem with = 0o~ -. Denote by e — 1/n, we 

can rewrite the Hill function as 

As a technical preparation, we now describe some elementary properties of the above 
specific Hill function. 

Lemma 4.1. Assume that e is sufficiently small. We have 

(a)Ify> I -J , then 

0{y) < Q e, y(3(y) < o e 

and if < y < e e , then 



0o > P(y) > - e) , and \y/3(y) - (3 y\ < /3b e. 



(4.2) 



(b )Moreover, 



d(yp(y)) 



and 



dy 

d(y0(y) - o y) 



< o e , for y > 



2e 



dy 



< o e , for < y < (j^) 6 ' 



Proof (a). If y > - , then 



e/(l-e) 



at \ 00 ^ 00 ^ 00 ^ Q 

m = TTT^ < 7^ < / n i/d- S ) < 



and 



v0M = ttt^ < ^ < Po£ - 



If < y < e £ , then 



> = TT^TTe > y 1 ' 6 ) ^ - e), 

1 + y L i £ 
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and 



1 + y L ' e 



(b) If y > (l/e) 2£ , then 



\(yp(y))'\=Po 
Since the function 



(1 - l)yV* 



(l+yl/e)2 



< /3 (- - l)y- 1,e < 



£. 



In 1 9 

(1 + -)X + -.T 2 

/(*) = " 



1+X 



e 2 , . £ 2 



is increasing for a; € (0, ) and f ( ) < e, then 

(1 + V e + V /S 
- /3 y) | = A) £ , , 1/e £ < fa, 

i{o< y< (-f-y.D 

1 + £ 

Returning to equation H4.1J1 . we let the initial function <p(t) be chosen in the 
following closed convex sets 

A(rj) = {(p(t) e C([0, 1]) : 1 + r? < <p(t) and (p(0) = 1 + 77}, 

where 77 is a small positive constant defined in the previous section. For given ip{t) in 
A(rf), we can have a unique solution to equation 114.111 . The relations 

Fp(t, <p) = y t , yt = y(t + s), -r < s < 0, t > 

define a continuous semiflow F = Fp on C([— r, 0]). 

We find that for the simpler equation 13.10 . if £ ^4(?7), then the solution 
will return to A(rj) after finite time. We like to know whether or not this situation 
still happens for equation 14, 1|) . The study of this point becomes necessary and also 
important in order to build a returning map. Fortunately, we have 



Lemma 4.2. Let y(t) be the solution of equation U-l\ ) with any initial function 
tp e A(rf). Then 

y{t)=x(t)+0(e\oge), 

fort £ [0,T X ] where T x is the period of periodic solution x(t), defined in Theorem \H.l\ 
to equation 
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Proof From J4.H we know that the solution y(t) is decreasing in t in the right 
neighborhood of the starting point t = 0. We can further claim that there exist three 
points r/i, t\, ij 2 , r\\ < if < 772, so that 

vim) = (J) 2 " > y(tl) = i, vim) = (^) , (4.3) 

v 2s 



and in the interval (0,772), the solution y(t) is decreasing. Indeed, if y(t) > I -J > 
(i) £/(1 ~ £) , and y(t - r) > ( - J > (I) £/(1 ~ e) , then by LemmaEHwe have 

< A>£, P{y T )y T < Po? 
and it follows from equation Il4.ll! that 

^L = -(S + P(y))y + kf3(y T )y T , 

< — - as e -> 0, (4.4) 

which means that j/(t) is decreasing and there exists a point 771 such that y(r?i) = 
(1/e) and I + 17 > > (1/e) , for i € (0,771). Similarly at the right neighborhood 
of 771, say (771, 771 + r/2),we have (3(y T )y T — 0(e) and ^ 14.41) still holds. This means the 
solution is still decreasing in t and there exist two points t\, r\%, r\\ < if < 772 so that 



y(*i) = i) y(m) 



e 2 



l+e 

By the Mean- Value Theorem, it is easy to know that 



or equivalently 



Therefore, 



\y(m) - y{m)\ >-x\m-m\ 



2 

m-m< -rivim) - vim)) = O(-eioge). 



t\-r}x < T72 — 771 = 0(-e log e). (4.5) 
Now using again the equations and \?>.il for t S [0, 771], we have from Lemma l4~D 

(x-y)' = -5(x-y)+0(e), 

which implies that 

\x(t) - y(t)\ = 0(e) 
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for t e [0, 771]. In particular at the point t = r]i, 

x(m) = Vim) + 0(e) = - + 0(e) = 1 + O(-eloge). 



It follows from l|3.2jl and l|3.4jl that t\, denned in (13 .3|) satisfies 

ti = 7/i + O(-eloge), 

and 

h =t\ + 0{-e\oge). 

For t € [7/1,7/2], since both the derivative of x(t) and y(i) are of the order of 0(1) 
and the length of the interval [771,772] is of the order of 0(— eloge), we can conclude 
that 

y(t)=x(t) + 0(-elog£). (4.6) 

For t e [772, t + 771], y(t -t)> (1/e) . By Lemma IQ1 we still have 

y(t-r)P(y(t-T)) = 0(e). 

By equation H4.1J1 we know that the solution y(t) will still be decreasing for t E 
[7/2, t + 7/1]. Note that < y(t) < 2/(772) = £ e , so that 114.211 in Lemma l4~D holds. Thus 
we can derive from (14. l|i that 

y'(t)=-ay(t) + 0(e), (4.7) 

for t 6 [772,^ + 771]. Coupling this equation with 1)3.411 and using 114.611 at the point 

t = 772 gives 

y{t)=x{t) + 0{e\oge) 

for t E [772, t + 7/1]. 

For i 6 [t + 771, r + 7/2], using again the fact that both the derivatives of x(t) and 
t/(£) are bounded by O(l) and the length of this interval is of order 0(— eloge), we 
have 

y(t)=x(t)+0(e\oge). 

For t > I + 7/2, the solution y{t) begins to increase since T satisfies l|3.7fl . By the 
similar argument used above, it follows that there exist three point 7/3, t\, 774, 7/3 < t\ < 
7/4 such that 

and 

7/3 =4 + 0(elog£), 7/ 4 = ^ + 0(eloge), (4.S 
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and 

t v 2 =t 2 + 0{e\oge). (4.9) 

We can continue the process above and it is not difficult to find that y{t) will satisfy 

y{t) =x(t)+0(eloge), 

for t £ [0,774]. 

By 113. 7J1 and 14.lt we find that the solution y{t) is increasing at the point t = r/4 
and in the interval [774, t + 773], y(t) and x(t) satisfy 

{x-y)' = -S(x-y)+0(e), 

and 

{x-y)\vi = O(eloge). 

So we have 

x-y = 0(e log e), for t g [774, r + 773] . 

For t £ [r + 773 , r + 774] , using the same argument as in the interval [r + 771 , r + 772] , we 
have 

y(t)=x(t)+0(e\oge). (4.10) 

Finally for t > t + 774, the solution is decreasing and attains the value 1 + 77 at some 
point T y . To be more specific, we have 

y(t)' = -6y(t) + 0(e), (4.11) 

and 

x{t) = -5x(t). (4.12) 
Using H4.10H . H4.11jl and Q4.12H we can derive that 

y(t) =x(t)+0(e\oge), (4.13) 

and 

T y =T x + 0{s logs). (4.14) 

Furthermore, we also have y{T y ) = 1 + 77 and 

y(t) > 1 + 77 for [T y - r, T y ] (4.15) 

provided that e is sufficiently small.D 

Remark 4.3. By Lemma \4-2\ and equation U-l\) we can have two positive con- 
stants Mi and M2 which are independent of e, so that 

\y(t)\<M u (4.16) 
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I < M, (4.17) 



Now we are ready to define a continuous returning map 

R : A(rf) 3 <p ^ y q{v) = Fp(q(<p),<p) e A(rj), 

where q(tp) — T y . In order to verify that there exists a unique fixed point in A(r)) for 
this map R, we need to derive its Lipschitz constant estimation and show this map R 
is contractive, i.e., its Lipschitz constant is less than 1. 

4.2. Lipschitz constant for the map R. Lipschitz constants of maps T : 
Dt — > Y, Dt Cl,l and Y normed linear space, are given by 

xeD T ,yeD T ,x^y \\X - V\\ 

In the case when Dt = X — R and a — [xi, X2] £ R, and / = T we set 

£[*i,z 2 ](/) = £(/|[zi,Z2]). 
In the case when / = yfi{y) 1 we define the following four Lipschitz constants 

L 2 = L 1 (y(3(y)), 

[(-) 2 *,+oo) 



£3= £(o,+oo)(y/%))) 



1 + e 



Similarly for the function / = y/3(y) — floy, we also define the following Lipschitz 
constant for later use, 

L 5 = L / f p2 \ =\ {y/3(y) — Poy)- 



l + e 



Theorem 4.4. When e is small, the Lipschitz constant for the map R satisfies 

lim Lr = < 1 

£^0 
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Proof Step 1. For 4>, 4> in A(rj). By a similar manner as in the proof of Lemma 
14.21 we conclude that there exist 771,772 and 771, 772 such that, respectively, 



V 2s 



and 

y ^ l) = (i) ' ^ fe) = ( TT7 )£ ' m-V2 = 0(-elag E ). 
Let 

Vmin = min{r?i,^i}, 

and 

?7max = max{?72,r?2}- 

Then by 114. 5 jl we have 

Vmax - Vmin = O(eloge). (4.18) 

For t £ [0, rfciin], using the equation l14.lt for y^(t) and , respectively, gives 

dyHt ) 
dt 

and 



= -[S + p(y+(t))]v*(t) + WivHt - T))y+(t - r), (4.19) 



= -[£ + + fc/%*(t - r))i/(t - r). (4.20) 

Now we begin to estimate the difference between y^(t) and y^(t). Coupling with l|4.19fl 
and H4.20H yields 

= 2/*) (4-21) 
-|/%V -/?(/)/] 
+fc[/3(^)^-/3(yf)j/f]. 

Substituting the following inequalities 

\P(v*)v+ - P{v*)v*\ <L 2 \y*-y% 

and 

\P{yt)yt - P{vf)vi\ < ^110 -0|| 

into (|4 . 2 If) . we have 

(/-/)' < (^ + L 2 )|/-/|+fcL 1 ||0-^||. (4.22) 
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Integrating H4.22J1 from to t, gives 

f ((5 + £ 2 )|^-/| + ^ip-^|)ds. 



ds. 



'0 

Similarly we have 

-(»* - V*) < | ((<* + £2) 1/ - y+\ + kL^ - 0||) 
Thus, 

1/ - /I < f ((* + L a ) 1/ - V*| + - 0ll) ds. 

Solving 114.2. ill ( or by Gronwall inequality), we obtain 

\y*-v*\ < Ci\\<l> - <f>\\, 

where 

e (5+L 2 )r, mi „ _ 1 

Ci = — feLx. 

+ L 2 

Step 2. Next for t £ [r] m m, ^max], we have 

I <L 3 |/-/|, 

and 

\P{yt)yt - P{yi)yt\ <Ii\\4>-4>\\- 
Thus from Q4.19I1 and 14.2011 we can obtain as before 

</ ((5 + L3)|^-/| + fcL 1 ||^-0||)d S + C 1 ||^- 

Then by Gronwall inequality we have 

I/- /I < cr a ||0- £|| 

where 

„5+L3(»?max-J)min) _ 1 

C 2 = Cie {5+L3){r '—- T '^ ) + - — kin > Ci. 

o + L 3 

Step 3. For t G [t? 

max; "7" T" Tjxxiin\i 

\0(v+)v+ - /3o/ - {P(V+)V* ~ PoV*)\ < L 5 |/ - y*\ 

and 

\0(yt)yt - P(yt)yt\ < L 2 \\<j> - 4>\\- 
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It is thus easy to have 

1/ - /I < £ ((<* + L s )\y+ - y*\ + kLzCrU - ds + C 2 ||«£ - 4>\\, 

and to conclude that (due to 1 + ri m i n — ?7max < t) 

|/-/|<C 3 ||0-^||, (4.28) 

where 

olt+tL 5 _ 1 

C 3 = C 2 e ar+TL5 + kL 2 d > C 2 . (4.29) 

a + L 5 

Step 4. When t > 1 + r/ m i n , We note that from the proof of Lemma [4.21 there 
exist 773 < 774, and 773 < so that 

= (j^)*' V^i) = ,VA~Vs = O(-eloge) 

and 

/fa) = (^Y, y*(m) = 0) , - 773 = o(- E iog e ). 

Let 

7 /min = mm{r7 3 ,77 3 }, ^ ax = max{?74, 774}. 
Then by 114. 8 jl we have 

Vmax-Vim = 0(elog £ ). (4.30) 

For t 6 [r + ?7mm,?7mm]i we similarly have 

l/"/l< / ((« + i5)|/-/|+fcL3C73||0-0||)d S + C 3 ||0-0||, 

and 

<C 4 ||0-^||, (4.31) 

where 

r w , , p( a + i 5)(')Ln- T -'7min) — 1 

C 4 = Cse^+^^-r-n^) + fcL3C3 > (4 _32) 

a - L 5 

Step 5. For t e [VmmiVtaxl 

\V+-V*\< /* ((5 + L 3 )\y*-y^\ + kL 4 C 4 \\^~4>\\)d S + C 4 \\<j>-^\\. 
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Thus, 

\V+-V*\ <C 5 \\<i>-$\\ (4.33) 

where 

e s+Ls(v^-vLj + £ i fci4 (4.34) 

<5 + ^3 y 

Step 6. For t € [7^ lax , T + lij. we have 

\V+ ~ V*\ < J ((<* + L 2 )\y* - y*\ + kL 3 C 5 \\<P - $\\) ds + C 6 \\<f> - $\\ 

and 

\y*-V*\<C e \\tt>-^\\, (4.35) 

where 

JS+L 2 )t _ i 

C 6 ^C 5 (e^ +L ^+ kL 3 ). (4.36) 

— L*2 

Step 7. When t > r + rj^ax) both y and y are decreasing and will take the value 
1 + r\ after finite time. Suppose that s and s satisfy 

2/ O) = 1 + V> = 1 + 17- 

For the later proof, we only consider the case s < s, since the case when s > s can be 
similarly dealt with and the proof will be omitted. By Lemma l4~2l we know 

y+(t) = x{t) + 0(e\oge), y*(t) = x{t) + O(eloge). 

By H4.8II . 114.911 and <I4 . 1 4|> we can also obtain 

s - (r + i£„) - T x - (r + i 2 ) + 0(e log e) 

and 

s-(t + r7n iax ) = T x - (t + t 2 ) + 0(e log e) 

where T x is the period of function x(t). Because the distant between r+T^ax and s may 
be greater than r. Thus we need split the interval [r + ^ al , s] as [T + r/^ax) 2T + f74ax]> 
[2r + ??4ax' 3r + J7^ ax ], • ■ ■ , [wT + 7^ ax , s] where the length of each interval is exactly r 
except the last one. Here m is the largest integer less than or equal to (s— (T+^ax))/ 7 "' 
We can estimate — &r| interval by interval and finally to obtain 

|l/*-V*|<CV||^-0||, ( 4 - 37 ) 




with 



(4.38) 
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For t G [a, a], the function y^(t) satisfies 

y*(t) = l + r] + O(eloge), t G [s, s] and y^(s) = 1 + rj, 
from which and the equation H4.1fl we know that y^(t) is decreasing and 

l^^l = \-(8 + P{y*))y* + kl3(yt)v} 

= \-6(l + v ) + 0(eloge)\, 
6(1 + Tj) 
2 

when e is small. Applying the Mean- Value theorem to the function y^(t) implies the 
existence of p G [a, a] such that 



|^(a) -</(a)| = |^(p)'(a-a)| > 



or 



by 



|s - s\ < g£^\v*W - = soT7) ly4>{s) ~ y " (s)l ' (4 - 39) 

Our ultimate goal is give the estimate of \yf (8) — yf(0)\ where 8 G [— r, 0]. Indeed 



\yl(6) ~ yt(0)\ < \vtifi) - yf(8)\ + \yf(0) - yf(8)\. (4.40) 
The first term of the right hand side is bounded by 

' s+e dy*(t) 



s+e dt 



-dt < M 2 a-a 



where Mq. is the maximum value of derivative of the function see remark 4.3. 

The second term of (14.4011 is bounded by CV||</> — (f>\\. Thus from H4.40I) . we have 

\yt(d) - yf(8)\ < C 7 (l + J^fiM ~ 4>\\- (4-41) 

When e — > 0, we have 

Li = 0( - ± U/ J , L 2 = 0(e), L 3 = 0(l/e), = 0(1), L 5 = 0(e). 
e(l + r;) 1 /^ 

Since L\ is exponentially small as e — > 0, we conclude from H4.25JI . 114, 27|) . i|4,29|l . 
(I4.32j) , ( I4.34j) , l|4.:ffijl and l|4.38jl that L# is exponentially small and satisfies 

lim L fl = lim C 7 (l + r/1 2 J = 0. 

£^0 £^0 V 6(1+7])' 

This completes our proof.Q 

Since Lr < 1, it means that the returning map R is contractive and there exists 
a unique fixed-point 4> in A(rj) so that i?(</>) = <f>- Thus we find a slowly oscillating 
periodic solution y(t,4>) for the equation 14. H . This periodic solution is also globally 
attractive for any initial function ip in A(rf). 
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5. Full Asymptotic Expansion for the Periodic Solution. In the previous 
section we use fixed-point theory to prove that there exists a unique periodic solution 
in ^4(77) for the equation y.lll . Now we like to give a quantitative analysis to this 
solution. Since the map R is contractive with the Lipschitz constant Fr being expo- 
nentially small, it is possible to give a full asymptotic expansion for this particular 
solution. 

If we take the initial function (j) = 1 + 77, then we get a solution y(t, 1 + 77) which 
may not be periodic. But by Lemma 4.2, y(t, l + rj) = x(t) + 0(e loge) when t is finite 
and there exists a T 1+v > such that 

y(T 1+ri , 1 + 77) = 1 + 77, y(T 1+v + 9, 1 + 77) > 1 + 77, 6 e [-r, 0). 

Then y(6 + T 1+v , 1 + 77) € A( V ). 

Assume that y{t) is the periodic solution to i|4.1Jl which can be extended to 
(—00, +00) and satisfies y(t) G A(rj) for t e [— r, 0]. Suppose also that for t G [— r, 0], 
y{t) has the following asymptotic formula 

00 

»(*)=!>(*) (5-1) 

i=0 

where <j>o{t) = y(t + Ti +7) ,1 + 77) and 4>i(t,e),i > 1 will be determined later. Let 
Vo{To + 0) denote the image of the returning map of R on (j>o, i.e., 

y o {T o + 6)=F {T o ,fo), 8e [-r,0] 

where T > satisfies 

yo(T ) = 1 + 7?, y (T + 9) > 1 + 77, 6 e [-r, 0). 
similarly by induction set 

MO) = yo(T + 6)-^ 

i-1 

yi-i(T i - 1 + e) = Fp(Ti-i,^2<f> j )=F fi (T i - 1 ,yi- 2 ), 1 > 2, 

and 

<k(0)=Vi-i(P + Ti-i)-yi-2{6 + Ti-2), 1 > 2, 
where for i > 1, Tj > and also satisfies 

W (Ti) = 1 + 77, w (Ti + 0) > 1 + 77, 6 G [-r, 0). 
Using the result in H4.39H we can also have 

\Ti-TiM < (Lny-'m-Tol 
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which means that the series 

oo 

is absolutely convergent to some constant, say T £ . Since by Theorem 14.41 we have for 
0€[-t,O], 

+ Ti) - yi _ x {6 + Ti_ x )\ = \F(T h w _ x ) - F(Ti_i, j/i_2)| (5.2) 

< L fi |y 4 _i(0 + I*-!) - W _ 2 (fl + T_ 2 )| 

< (^) l - 1 |j/i(^ + ri)-j/o(e + T )| 

< (iij)iy o (t,0o)-^o|. 

Then for s <E [— r, 0], we conclude that 

oo 

y (s + T ) + ]T (vA s + t j) ~ Vj-i(* + T j-i)) (5-3) 

3=1 

is absolutely convergent and have the same value as the initial function S^lg </>i(s). 
Since Lr is exponentially small, It is easy to prove that the value of T £ is dominated 
by T and likewise the value of l|5.3J) is dominated by yo(s + T ), each of which has the 
exponential error bound. For the leading term j/o(t), obviously we have from Lemma 
14.21 the rough estimates 

y (t) - x{t) + O(eloge), T = T X + 0(e\oge). 

Next we would like to give the refined estimate for yo[t) and To by using this infor- 
mation and developing the idea of Lemma l4~2l 

Like in Lemma l4~2l we shall split the interval [0,T ] and estimate yo(t) interval 
by interval. For illustration, we only need to show the first interval's estimate to the 
readers. Keep in mind that the initial data is taken as y(t + T x+V , 1 + Tj) which is 
greater than 1 + r\ when t lies in the interval [— r, 0). Let 771 and 772 be still the value 
in Lemma 4.2 and t\° satisfy yo(tf°) — 1. We denote yo{t) = yo(t). Integrating the 
equation 14.lt from to t, t £ [0, if ], gives 

Vo(t) -W)(Q) = -S f y (t)dt- [ (3(y )y (t)dt (5.4) 
Jo Jo 

+k f (3(y (t-T))y (t-T)dt. 
Jo 

It is obvious that the last term of the right hand side of llo.4ll is exponentially small 
and can be viewed as O(e). Next we claim that 



/ P(y )y (t)dt = O(e). 
Jo 



(5.5) 
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To see this, we need to note that yo(t) is decreasing in t € [0,if] and dy /dt is the 
order of O(l) or more precisely 

-a(l + V) < ^ = -\P(Vo) + %o + 0(e) < -6 + O(e). (5.6) 
Thus from H5.6J1 and the fact 

t y ° 

P(y )yo(t)^dt < I" (3(y )y (t)^dt, 



dt J n dt 



f 1 p y 



1 1+7} 

= 0(e), 

we know that J (3(y )y (t)dt is also the order of 0(e) and the claim of l|5.5|) is true. 
It follows then from l|5.4jl 



y (t) =-S y (t)dt+l + r) + O(e). 
Jo 

Solving this integral equation, we get 

y Q (t) = (l + V + 0(e))e- St , 

which implies 

yo(t) = y (t, e) = x{t) + 0(e), t e [0, t\]. (5.7) 

Using the same approaches above we can prove that in the whole interval [0, To], Q5.7J1 
is still true. Furthermore, we can prove that 

To = T x + O (e) , 

which completes our refined estimates. 



REFERENCES 



[1] M. Adimy and L. Pujo-Menjouet, A mathematical model describing cellular division with 

a proliferating phase duration depending on the maturity of cells, Elect. Journal of Din". 

Eq., 107 (2003), pp. 1-14. 
[2] J. Belair, M. C. Mackey, and J. M. Mahaffy, Age- structured and two-delay models for 

erythropoiesis, Math. Biosci., 128 (1995), pp. 317-346. 
[3] S. Bernard, J. Belair, and M. C. Mackey, Sufficient conditions for stability of linear 

differential equations with distributed delay, Discrete Contin. Dyn. Syst. Ser. B, 1 (2001), 

pp. 233-256. 

[4] A.-M. Buckle, R. Mottram, A. Pierce, R. N. Lucas, G. S., J. A. Miyan, and A. D. 

Whetton, The effects of bcr-abl protein tyrosine kinase on maturation and proliferation 

of primitive haemotopoietic cells, Molecular Medicine, 6 (2000), pp. 892-902. 
[5] F. Burns and I. Tannock, On the existence of a Go phase in the cell cycle, Cell Tissue 

Kinet., 3 (1970), pp. 321-334. 
[6] C. J. Eaves and A. C. Eaves, Stem cell kinetics, Baillieres Clinical Haematology, 10 (1997), 

pp. 233-257. 

[7] P. Fortin and M. C. Mackey, Periodic chronic myelogenous leukemia: Spectral analysis of 
blood cell counts and etiological implications, Brit. J. Haematol., 104 (1999), pp. 336-345. 



M. C. MACKEY, C. OU, L. PUJO-MENJOUET and J. WU 



T. Haferlach, M. Winkemann, C. Nickenig, M. Meeder, L. Rammpetersen, 
R. Schoch, M. Nickelsen, K. Webermatthiesen, B. Schlegelberger, C. Schoch, 
W. Gassman, and H. Loppler, Which compartment are involved in philadelphia- 
chromosome positive chronic myeloid leukaemia? an answer at the single cell level by 
combining may-grunwald-giemsa staining and fluorescence in situ hybridization techniques, 
Brit. J. Haematol., 97 (1997), pp. 99-106. 

C. Haurie, D. C. Dale, and M. C. Mackey, Cyclical neutropenia and other periodic hema- 
tological diseases: A review of mechanisms and mathematical models, Blood, 92 (1998), 
pp. 2629-2640. 

T. Hearn, C. Haurie, and M. C. Mackey, Cyclical neutropenia and the peripheral control 
of white blood cell production, J. Theor. Biol., 192 (1998), pp. 167-181. 

X. Jiang, Eaves, C. J., and A. C. Eaves, IL-3 and G-C'SF gene expression in primitive 
PH+C'D34(+) cells from patients with chronic myeloid leukemia ( CML), Blood, 90 (1997), 
pp. 1745-1745. 

L. Kold-Andersen and M. C. Mackey, Resonance in periodic chemotheray: A case study 
of acute myelogenous leukemia, J. Theor. Biol., 209 (2001), pp. 113-130. 

L. G. Lajtha, On DNA labeling in the study of the dynamics of bone marrow cell populations, 
Grune & Stratton, New York, 1959, pp. 173-182. 

M. C. Mackey, Unified hypothesis of the origin of aplastic anemia and periodic hematopoiesis, 
Blood, 51 (1978), pp. 941-956. 

, Dynamic haematological disorders of stem cell origin, in Biophysical and Biochemical 

Information Transfer in Recognition, J. G. Vassileva-Popova and E. V. Jensen, eds., Plenum 
Publishing Corp., New York, 1979, pp. 373-409. 

, Periodic auto-immune hemolytic anemia: An induced dynamical disease, B. Math. 

Biol., 41 (1979), pp. 829-834. 

, Mathematical models of hematopoietic cell replication and control, in The Art of Math- 
ematical Modeling: Case Studies in Ecology, Physiology & Biofluids, H. G. Othmer, F. R. 
Adler, M. A. Lewis, and J. C. Dalton, eds., Prentice Hall, 1997, pp. 149-178. 

M. C. Mackey, Cell kinetic status of hematopoietic stem cells, Cell Prolif., 34 (2001), pp. 71- 
83. 

M. C. Mackey and P. DSrmer, Continuous maturation of prolif erating erythroid precursors, 

Cell Tissue Kinet., 15 (1982), pp. 381-392. 
M. C. Mackey and R. Rudnicki, Global stability in a delayed partial differential equation 

describing cellular replication, J. Math. Biol., 33 (1994), pp. 89-109. 
C. Ou and J. Wu, Periodic solutions of delay differential equations with a small parameter: 

Existence, stability and asymptotic expansion, (2004). 
L. Pujo-Menjouet, S. Bernard, and M. C. Mackey, Long period oscillations in a go 

model of hematopoietic stem cells. 
L. Pujo-Menjouet and M. C. Mackey, Contribution to the study of periodic chronic myel- 
ogenous leukemia, (2004). 
S. I. Rubinow and J. L. Lebowitz, A mathematical model of neutrophil production and 

control in normal man, J. Math. Biol., 1 (1975), pp. 187-225. 
M. Santillan, J. Belair, J. M. Mahaffy, and M. C. Mackey, Regulation of platelet 

production: The normal response to perturbation and cyclical platelet disease, J. Theor. 

Biol., 206 (2000), pp. 585-903. 
J. A. Smith and L. Martin, Do cell cycle?, Proc. Natl. Acad. Sci. USA, 70 (1973), pp. 1263- 

1267. 

J. Swinburne and M. C. Mackey, Cyclical thrombocytopenia: Characterization by spectral 
analysis and a review, J. theor. Med., 2 (2000), pp. 81-91. 

N. Takahashi, I. Miura, K. Saitoh, and A. B. Miura, Lineage involvment of stem cells 
bearing the Philadelphia chromosome in chronic myeloid leukemia in the chronic phase 
as shown by a combination of fluorescence-activated cell sorting and fluorescence in situ 
hybridization, Blood, 92 (1998), pp. 4758-4763. 

H. O. Walther, Contracting return maps for monotone delayed feedback, Discrete and Con- 
tinuous Dynamical Systems, 7 (2001), pp. 259-274. 

, Contracting return maps for some delay differential equations, Topics in Functional 

Differential and Difference Equations, Fields Institute Communications, (T. Faria and P. 
Freitas eds.) Providence: American Mathematical Society, 29 (2001), pp. 349-360. 

, Stable periodic motion for a system with state dependent delay, Differential and Integral 

Equations, 15 (2002), pp. 923-944. 

, Stable periodic motion of a system using echo for position control, preprint, (2003). 

G. F. Webb, Theory of Nonlinear age-dependent population dynamics, vol. 89 of Monographs 



PERIODIC SOLUTIONS IN CHRONIC MYELOGENOUS LEUKEMIA 



23 



and Textbooks in Pure and Appl. Math. Series, Dekker, 1985. 
[34] J. Wu, Stable phase-locked periodic solutions in a delay differential system, 3. Differential 
Equations, 194 (2003), pp. 237-286. 



